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Abstract. In this paper we present an algorithmic approach to the generation of fully con- 
servative difference schemes for linear partial differential equations. The approach is based 
on enlargement of the equations in their integral conservation law form by extra integral 
relations between unknown functions and their derivatives, and on discretization of the ob- 
tained system. The structure of the discrete system depends on numerical approximation 
methods for the integrals occurring in the enlarged system. As a result of the discretization, 
a system of linear polynomial difference equations is derived for the unknown functions and 
their partial derivatives. A difference scheme is constructed by elimination of all the partial 
derivatives. The elimination can be achieved by selecting a proper elimination ranking and 
by computing a Grobner basis of the linear difference ideal generated by the polynomials 
in the discrete system. For these purposes we use the difference form of Janet-like Grobner 
bases and their implementation in Maple. As illustration of the described methods and al- 
gorithms, we construct a number of difference schemes for Burgers and Falkowich-Karman 
equations and discuss their numerical properties. 
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1 Introduction 

It is well-known that finite differences along with finite elements and finite volumes are most 
important discretization schemes for numerical solving of partial differential equations (PDEs) 
(see, for example, H El H El El 13 El)- 

Mathematical operations used in the construction of difference schemes for PDEs are sub- 
stantially symbolic. Thereby, it is a challenge for computer algebra to provide an algorithmic 
tool for automatization of the difference schemes constructing as well as for the investigation of 
properties of the difference schemes. One of the most fundamental requirements for a difference 
scheme is its stability which can be analyzed with the use of computer algebra methods and 
software 0. 

Furthermore, if PDEs admit a conservation law form or/and have some symmetries, it is 
worthwhile to preserve these features at the level of difference schemes too. In particular, 
a tool for automatic construction of difference schemes should produce conservative schemes 
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whenever the original PDEs can be written in the integral conservation law form. One of 
such tools GRIDOP written in Reduce jlUl 111) is based on symbolic operator methods and 
generates conservative finite-difference schemes on rectangular domains in an arbitrary number 
of independent variables. However, the generation is not entirely automatic. A user of GRIDOP 
has to specify function spaces together with associated scalar products and define grid operators 
as finite-difference schemes. Then the user may provide partial differential equations in terms 
of the defined grid operators or the adjoints of those operators. Under these conditions the 
package returns the finite-difference equations for the dependent variables. 

Besides, a few other applications of computer algebra are known to construct finite-difference 
schemes |12U13| which, being also not completely automatic, are applicable to PDEs of a certain 
form. 

In this paper we describe a universal algorithmic approach to the automatic generation of 
conservative difference schemes for linear PDEs with two independent variables admitting the 
conservation law form. This approach generalizes and extends the observations of paper |14j 
where it was noticed that a conservative difference scheme can be derived as a compatibility 
condition for a system of difference equations. The system is composed of a discrete form 
of the original PDEs taken in the integral conservation law form and of a number of natural 
integral relations between functions and their partial derivatives. The finite-difference scheme 
is obtained by elimination of all the partial derivatives from the system. We also show, by the 
example of Burgers equation, that one can also apply the difference elimination approach to 
generate of difference schemes without use of conservation law form. 

To perform the difference elimination we apply the Grobner bases method invented 40 years 
ago by Buchberger |15| for polynomial ideals. This method has become the most universal 
algorithmic tool in commutative algebra and algebraic geometry and found also numerous fruitful 
applications for computations in certain noncommutative polynomial rings as well as in rings 
of linear differential operators and differential polynomials [TJI|- Nowadays, all modern general- 
purpose computer algebra systems, for example, Maple Jlj and Mathematica jTHj, have special 
built-in modules implementing algorithms for computing Grobner bases. However, the fastest 
implementation of these algorithms for commutative polynomial algebra is done in the special- 
purpose systems Singular ^5] and Magma [20]. As to the difference algebra |U, in spite of 
known for long time (see [22] and references therein) extensions of Buchberger's algorithm [22] to 
difference polynomial rings, there are only a few implementations of the algorithm specialized to 
shift Ore algebra: in the Ore algebra library package of Maple [23] , in the library OreModules [2E] 
developed using the latter package and in Singular (Plural) j2B] . These packages can be used for 
computing Grobner bases of linear difference ideals and modules, and, in particular, for those 
linear systems which are considered below. 

In the given paper we present, however, another algorithm for computing difference Grobner 
bases. This algorithm is superior over our Janet division algorithm whose polynomial version |27| 
in most cases is computationally more efficient than Buchberger's algorithm |2B]- In addition, 
unlike the above mentioned implementations of Buchberger's algorithm for the shift Ore algebra, 
the algorithm described below and its recent implementation in Maple [22] admit a natural 
extension to nonlinear difference systems exactly in the same way as differential involutive 
algorithms |30[ I31| I32j . The algorithm improves our Janet-like division algorithm [SS] adapted 
to linear difference ideals [33] • The improvement includes, in particular, the difference form of 
the involutive criteria [2Z] modified for Janet-like reductions. These criteria allow to avoid some 
useless reductions, and thereby accelerate the computation. 

The structure of the paper is as follows. In Section 2 we describe the basic idea of our 
approach to the generation of finite-difference schemes for two-dimensional PDEs. Section 3 
contains definitions and notions of difference algebra which are used in the sequel. Here we 
define Grobner bases for linear difference ideals and their special form called Janet-like bases. 
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In Section 4 we present an improved version of the algorithm in paper [31] and briefly discuss 
some relevant computational aspects. Section 5 illustrates our approach to construction of 
difference schemes by simplest second-order equations - Laplace's equation, the wave equation, 
the heat equation, and by the first-order advection equation. In Section 6 we generate several 
difference schemes for Burgers equation. But for all that, to avoid problems arising in computing 
of nonlinear Grobner bases, we denote the square of the dependent variable by an extra function. 
Besides, we characterize some of the constructed schemes by the modified equation method. In 
Section 7 we consider the two-dimensional quadratically nonlinear Falkowich-Karman equation 
describing transonic flow in gas dynamics. Here, we succeeded in computing of the nonlinear 
Grobner basis by hand, and in that way generated the cubic nonlinear difference scheme which 
possesses some attractive properties. These properties as well as those of the schemes generated 
for Burgers equation are illustrated by some numerical experiments in Section 8. We conclude 
in Section 9. 



2 Basic idea 

It is well-known [I] IE] that a rather wide class of scalar PDE and some systems of PDEs 
can be written in the conservation law form 

dv d „. , 

where v is a m-vector function in the unknown n-vector function u and its partial derivatives 
u x ,u y , u xx , .... The vector function F maps R m into R m . 

By Green's theorem (curl theorem in the plane), vector PDE Q is equivalent to the integral 
relation 



-F(v)dx + vdy = 0, (2) 

T 

where T is an arbitrary closed contour. Approximation of (J2J) rather than of Q on a difference 
grid (balance or integro-interpolation method) is a natural way to generate conservative finite- 
difference schemes for PDEs of order two and higher. 

Throughout this paper we shall consider orthogonal and uniform grids with the grid mesh 
steps hi and hi 

x j+1 - Xj = hi, y k+1 -y k = h 2 (3) 

and denote the grid values of the vector function u(x,y) and all its partial derivatives occurring 
in © by 

u(x,y) =^u(xj,y k ) = u jk , u x (x,y) =^u x (xj,y k ) = (u x ) jk , 

u y( x ,y) =^ u y(xj,yk) = (Uy)jk, u xx (x,y) u xx (xj,y k ) = (u xx ) jk , (4) 



Choose the integration contour as follows (see Fig.^) and add all the related integral relations 
between the dependent vector variable and its partial derivatives: 

xj+2 ryk+2 

u x dx = u(x j+2 ,y) -u(xj,y), / u y dy = u(x,y k+2 ) - u(x,y k ), 

r x j+2 ryk+2 

/ u xx dx = u x (x j+2 ,y) - u x (xj,y), / u xy dy = u x (x,y k+2 ) - u x (x,y k ), (5) 
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Figure 1. Integration contour on grid. 

To obtain a source system of discrete equations for constructing of a difference scheme, we 
consider numerical approximations of the integral equations (j2J for the contour of Fig. ^ and of 
the relations © in terms of the grid unknowns @. Although generally one can use different 
numerical approximations for the integral equations in (j2J and ©, we apply here for all these 
equations the simplest rectangle (midpoint) rule 



F{v) j+lk+2 



u 



x )j+l k 



F(v) j+lk + (v j+2k +i - v jk+ i) = 0, 
2hi=Uj + 2k~Ujk, ( u y)jk+i ■ 2^2 = ' 



j k+2 



U 



jk, 



(6) 



Thereby, we obtained system © of difference equations in grid unknowns (|4j). In doing so, 
the number of scalar equations added to the (vector) integral equation © corresponds to the 
number of proper partial derivatives of order less than or equal to the order of partial derivatives 
involved in the integrand of 

It follows that eliminating from © all the proper grid partial derivatives gives equations con- 
taining only independent (vector) function u, and, hence composing a finite-difference scheme. 
If system © is linear, then this difference elimination can always be algorithmically achieved 
by the Grobner bases method considered in the next section. 

It should be noted that generation of finite-difference schemes on grid by the elimination 
can be also applied to PDEs irrelative to their conservation law properties. Again, one has to 
add to the initial differential equations, written in terms of grid variables (|4"|). the corresponding 
number of integral relations © and approximate them by numerical quadrature formulas. Such 
an approach may give more flexibility in generation of distinct difference schemes, and we apply 
it in Section 4 to the first-order advection equation, and in Section 5 to Burgers equation. Gene- 
rally, however, for the secondhand higher-) order PDEs admitting the integral conservation law 
form, the difference scheme obtained directly from the differential form may not be conservative. 
Besides, the difference elimination based on the integral form is usually more efficient than 
that based on the differential form. This is because the number of partial derivatives to be 
eliminated in the former case is smaller than in the latter case. Indeed, the integrand in (J2J) has 
the differential order smaller by one than that in (^Q) whereas computational complexity of the 
elimination is at least exponential in the number of eliminated variables ■ 

3 Difference Grobner bases 

A difference ring R is a commutative ring with a unity together with a finite set of mutually 
commuting injective endomorphisms 8i,...,9 n of R. Similarly, one defines a difference field. 
Elements {y , . . . , y m } in a difference ring containing R are said to be difference indeterminates 
over R if the set 

• • • 0i V I {h,...,k n } e Z^ , 1 < j < m) 

is algebraically independent over R. 
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Hereafter we shall consider the ring of functions of n variables xi,...,x n with the basis 
endomorphisms B{ o f(x±, . . . , x n ) = f{x\, . . . ,Xi + 1, . . . , x n ). acting as shift operators. 

The field Q(xi, . . . , x n ) of rational functions in {x\, . . . , x n } whose coefficients are rational 
numbers is an example of difference field, and we shall assume in the next sections that the 
coefficients of PDEs belong to this field. 

Let K be a difference field, and R := IKjy 1 , . . . ,y m } be the difference ring of polynomials 
over K in variables {9^ o y k | /i E ^>o> k = 1, . . . ,m). Hereafter, we denote by Rl the set of 
linear polynomials in R and use the notations: 

n 

:= {91* | /x E Z|„}, de gi (^ o y k ) : = W , deg(0" o /) := := £ fi h 

i=i 

lcmO*,!/) :={max{ W ,^i},...,max{ Mn , !/„}}, lcm(r o y k , 9 V o y k ) := 9 lcm ^ o y k , 
0* o y k c 9 V o / when Z> A |i/ - /i| > 0. (7) 

A difference ideal is an ideal ICR closed under the action of any operator from 0. If F := 
• • • j fk} C R is a finite set, then the smallest difference ideal containing F will be denoted 
by Id(-F). If for an ideal I there is F C Rl such that / = Id(-F), then I is a linear difference 
ideal. 

A total ordering y on the set of 9^ o y 3 is a ranking if V i, j, k, fi, v the following hold: 
9 i 9 l * o yi y o y3, o y9 V oy k «=> 0*0^ o >- 0^" o /. 

If >- \u\ 9^ o yi y 9 U o y k the ranking is orderly. If j > k ==>■ 9^ o y J >- 9 U o y k the 
ranking is elimination. 

Given a ranking >-, a linear polynomial / E Rl \ {0} has the leading term of the form atfot/ 3 , 
i? E 0, where $oy J is maximal w.r.t. y among all 9 fl oy k which appear with nonzero coefficient 
in /. lc(/) := a E K \ {0} is the leading coefficient and lm(/) := i? o y J is the leading (head) 
monomial. 

A ranking acts in Rl as a monomial order. If F C \ {0}, then lm(F) will denote the set 
of the leading monomials and lmj(F) will denote its subset for indeterminate y 3 . Thus, 

lm(F) = \jp =1 ]mj(F). 

Given a nonzero linear difference ideal / = Id(G) and a ranking >-, the ideal generating set 
G = { gi , . . . , g s } C R L is a Grdbner basis [221 EI] of I if 

V/E/nR L \{0}, 3geG, 0E0 : lm(/) = 6 o lm( 5 ) . (8) 

It follows that the head monomial of / E I \ {0}, as well as the polynomial / itself, is reducible 
modulo G and yields the head reduction: 

f^f':=f-lc(f)9o(g/lc(g)), f' E I. 

9 

If /' ^ 0, then its leading monomial is again reducible modulo G, and, by repeating the reduction 
finitely many times |161 122[ 123] we obtain / — ► 0. Generally, if a polynomial h E Rl contains 

G 

a term with monomial u and coefficient c / such that u = •& o lm(/) for some i? E and 
/ £ F C Rl \ {0}, then h can be reduced: 

h^ti :=h-c9o(f/\ c (f)). (9) 

By applying the reduction finitely many times, one obtains a polynomial h which is either zero 
or such that all its (nonzero) terms have monomials irreducible modulo the set F C Rl. In both 
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cases h is said to be in the normal form modulo F (denotation: h = NF(h, F)). A Grobner 
basis G is reduced if 

VgeG : g = NF(g,G\{g}). 

In our algorithmic construction of reduced Grobner bases we shall use a restricted set of 
reductions called Janet-like (cf. [SHUSH) and defined as follows. 

For a finite set F C Ml and a ranking >-, we partition every ]mj~(F) into groups labeled 
by do, . . . , di G Z>o, (0 < i < n). Here [0]k := lnifc(F) and for i > the group [do, ■ ■ ■ , ck]k is 
defined as 

[d , . . . ,di] k := {u G lm fe (F) j d = 0,dj = deg^u), 1 < j < i}. 
Now we characterize a monomial u G lmk(F) by the nonnegative integer Aj: 

Xi(u,lm k (F)) := max{deg,(t») | u,v e [d , . . . ,di-i]fc} - degj(u). 
If Aj(n,lm fc (F)) > 0, then 0f such that 

Si := min{degj(t;) - deg^n) [ u,v G [d , . . . ,di-i]fc, degj(v) > degj(u)} 

is called a difference power for f G F with lm(/) = 

Let DP(f, F) denotes the set of all difference powers for / £ F. Now we define the partition 
of the set into two disjoint subsets 

Q(f,F):={e»\3e v eDP(f,F) : fi - v € Z^ }, J(f,F):=@\Q(f,F), 

which is similar to the partition of monomials into nonmultiplicative and multiplicative ones in 
the involutive approach |27j . 

A finite basis G of / = Id(G) is called Janet-like [3B1ISI] if 

V/£/nRi\{0}, 3geG,6e J(g,G) : lm(f) = 6o\m(g). (10) 

In full analogy with © a J -reduction is defined as 

h^h':=h-c9o(f/lc(f)), 9eJ(f,F), (11) 

for polynomial h G Ml containing monomial u with coefficient c ^ satisfying u = $ o lm(/) 
for some / € F C R L \ {0} and ■& G J(/, F). 

Apparently, any element in the ideal I = Id(G) is J -head reduced to zero by the finite 
sequence of jT-head reductions by elements g G G in the Janet-like basis G: 

f^f':=f-Hf)0o{g/]c{g)), 6 e J(g,G), f'eL (12) 

If the leading monomial of p G M \ {0} is not ^-reducible modulo a finite subset Fcl \ {0} 
we say that p is in the J- head normal form modulo F and write p = HNFj(p, F). If none of 
monomials in p is j7-reducible modulo F we say that p is in the (full) normal form modulo F 
and write p = NFj(p, F). 

Since j7-reducibility implies the Grobner reducibility Q, a Janet-like basis satisfying IjlOJ) is 
a Grobner basis. The converse is generally not true, that is, not every Grobner basis is Janet-like. 

The properties of a Janet-like basis are very similar to those of a Janet basis [SHI, but the 
former is generally more compact than the latter. For all that we consider hereafter only minimal 
bases. 
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Let GB be a reduced Grobner basis, satisfying [25] : 

VgeGB : g = NF(g, GB \ {<?}). (13) 

Let now J.B be a Janet basis, and J LB be a Janet-like basis of the same ideal and for the 
same ranking. Then for their cardinalities the inclusion Card(G-B) < Card( JLB) < Card(JS) 
holds [271 133] . Here Card abbreviates cardinality, that is, the number of elements. 

Whereas the algorithmic characterization of a Grobner basis is zero redundancy of all its 
5-polynomials |15l 123] . the algorithmic characterization of a Janet-like basis G has the following 
form (cf. [33]): 

VgeG, $£DP(g,G): NFj(ti o g, G) = 0. (14) 

These conditions are at the root of the algorithmic construction of Janet-like bases as described 
in the following section. 



4 Algorithm 

In this section we present an algorithm for constructing a reduced Grobner basis (|E]) of the ideal 
generated by an input set of linear difference polynomials. The algorithm is an improved version 
of the algorithm in paper [HI] and translates to the difference form of the polynomial involutive 
algorithm ]27] modified for the Janet-like reductions. 

To apply the difference form of criteria to avoid some unnecessary reductions we need the 
following definition. 

An ancestor of a difference polynomial / eFc Kl\{0} is a polynomial g £ F of the smallest 
deg(lm(<7)) among those satisfying / = 6 o g modulo Id(F \ {/}) with 9 E &■ If for all that 

deg(lm(5f)) < deg(lm(/)), 

then the ancestor g of / is called proper. 

If an intermediate polynomial h that arose in the course of the below algorithm has a proper 
ancestor g in the intermediate basis G, then h has been obtained from g via a sequence of shift 
operations of the form 'dog where $ € DP(g, G) with lm(i? o g) ^-irreducible modulo G. For 
the ancestor g itself the equality lm(anc(<7)) = lm(g) holds. 

In the main algorithm GrobnerBasis and its subalgorithms we endow every element / € G 
in the intermediate set of difference polynomials G (occuring in the set T) with a triple structure 
of the form: 

P = {/, 9, dpow} , 
where 

pol(p) := / is the polynomial / itself, 

anc(p) := g is an ancestor of / in G, 

dp(p) := dpow is a (possibly empty) subset of DP(f,G). 

The set dpow associated with the polynomial / accumulates all the difference powers for / 
which have been already applied to / in the course of the algorithm. Keeping this information 
serves to avoid useless repeated applications of the difference power operators. Knowledge of 
ancestors for elements in the intermediate basis helps to avoid some unnecessary reductions by 
applying Buchberger's chain criterion [23] adapted to Janet-like reductions. 
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Algorithm: GrobnerBasis (F, -<) 



Input: F G M.l \ {0}, a finite set; -<, a ranking 


Output: G, a reduced Grobner basis of Id(i ? ) 


1 


choose / G F with the lowest lm(/) w.r.t. >- 


2 


T:= {f,f,0} 


3 


Q:={{q,q,0}\qeF\{f}} 


4 


Q :=HeadReduce(Q, T, y) 


5 


while Q/0do 


6 


choose p Q with the lowest lm(pol(p)) w.r.t. >- 


7 


Q -=Q\ {p} 


8 


if pol(p) = anc(p) then 


9 


for all { q G T | lm(pol(g)) = 0^ o lm(pol(p)), > } do 


10 


Q:=QU{q}; T:=T\{q} 


11 


od 


12 


fi 


13 


h :=TailNormalForm(p, T, -<) 


14 


T := T U {/i, anc(p), dp(p)} 


15 


for all g G T and G £>P(g, T) \ dp(g) do 


16 


Q := Q U o pol(g), anc(g), 0}} 


17 


dp(g) :=dp(g)nL>P(g,r)U{i?} 


18 


od 


19 


Q :=HeadReduce(Q,T, -<) 


20 


od 


21 


return {pol(/) / G T} or {pol(/) | / G T | / = anc(/)} 



In the above main algorithm GrobnerBasis and its subalgorithms presented below, where 
no confusion can arise, we simply refer to the triple set T as the second argument in DP, NFj, 
and HNFj instead of the polynomial set {g = pol(t) | t G T}. Sometimes we also refer to the 
triple p instead of pol(p). Besides, when we speak of reduction of the triple set Q modulo triple 
set T we mean reduction of the polynomial set 



{/ = pol(g) | q G Q} 
modulo 

{g = pol(i) | t G T}. 

Correctness and termination of algorithm GrobnerBasis can be shown exactly as in the 
polynomial case |27l 133] . Here we only elucidate some related features of the algorithm. 

At steps 4 and 19 the jT-head reduction is performed for the difference polynomials in Q 
modulo those in T. Then the remaining tail reduction is done in line 13 to obtain the (full) 
^7-normal form. Thereby, the main while-loop 5-20 terminates when the conditions (|14|) hold 
for the difference polynomial set G composed from the first elements of triples in T 

G := {po%) | g G T}, (15) 

and the set Q is empty. The upper for- loop 9-11 provides minimality of the output Janet-like 
basis contained in T |27j . Another for-loop 15-18 constructs new conditions (|14|) which have to 
be further examined because of the insertion of a new element in T at step 14. Besides, the set 
of difference powers is upgraded at step 17. 
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Furthermore, the main algorithm GrobnerBasis together with its subalgorithms presented 
below ensures that every element in the output Janet-like basis composed from the first elements 
in the triple set T has one and only one ancestor. This ancestor is apparently irreducible, in the 
Grobner sense @, by other elements in the basis. Thereby, those elements in the output basis 
that have no proper ancestors constitute the reduced Grobner basis (|13|) that is returned by the 
main algorithm at the last step 21. 

The algorithm HeadReduce invoked in lines 4 and 19 of the main algorithm returns the 
set Q which, if nonempty, contains part of the intermediate basis ^7-head reduced modulo T. 
The reductions are performed by its subalgorithm HeadNormalForm that is invoked at step 6 
of the algorithm. 

If algorithm HeadNormalForm returns /i / 0, then \m{h) does not belong to the initial 
ideal generated by {lm(pol(/)) | / G Q U T} (23 EH]- In this case the triple {h, h, 0} for h is 
inserted (line 9) into the output set Q. Otherwise, the output set Q retains the triple p as it is 
in the input. 

In the case when h = and pol(p) has no proper ancestors that is checked at step 14, 
all the descendant triples for p, if any, are deleted from the intermediate set S at step 16. 
Such descendants cannot occur in T owing to the choice conditions at steps 1, 6 and to the 
displacement condition of step 9 in the main algorithm GrobnerBasis. Steps 14-18 serve for 
the memory saving and can be ignored if the memory restrictions are not very critical for a given 
problem. In this case all those descendants will be casted away by the criteria checked in the 
below algorithm HeadNormalForm. 

Algorithm HeadNormalForm performs verification (step 3) of jT-head reducibility of the 
input polynomial h modulo the polynomial set (|15|). This verification consists in searching 
a difference polynomial (reductor) in the set G defined in (|T5|) such that G yields the reduc- 
tion (|11|). If the search fails, that is, there is no ./-reductor, the algorithm returns at step 4 the 
input polynomial. 

For the ^7-head reducible input polynomial pol(p) that is checked at step 3 of algorithm 
HeadNormalForm, the following three criteria are verified at step 9 

Criteria(p,g) = d(p,g) V C 2 (p,g) V C 3 (p,g) , (16) 

where 

Ci(p,g) is true -4=^ lcm(lm(anc(p)), lm(anc(<?))) C lm(pol(p)), 

C-2 iPi d) is true 3 t € T such that 

lcm(lm(pol(i)), lm(anc(p))) C lcm(lm(anc(p)), lm(anc(g))) A 
lcm(lm(pol(t)), lm(anc((?))) C lcm(lm(anc(p)), lm(anc(p))), 

C 3 (p,g) is true <=^> 3 t € T A y € NM c (t,T) with lm(pol(i)) • y = lm(pol(p)), 
lcm(lm(anc(p)), lm(anc(i))) C lm(pol(p)) A idx(i, T) < idx(/, T), 
where idx(t, T) enumerates the position of triple t in set T. 

In aggregate, criteria (|16|) translate (cf. |27| I37j) Buchberger's chain criterion |23] into the linear 
difference algebra. 

In addition, if all difference polynomials in the input set F for the main algorithm Grobner- 
Basis have constant coefficients, then the set of criteria 1)16(1 can be enlarged with one more 
criterion C4: 

C^(p,g) is true for lm(pol(p)) = 6oy k , lm(pol(p)) = $ o y k <^=> \cm(9,'&) =9$. 
Criterion C4 is the difference form (cf. [21]) of Buchberger's co-prime criterion |23) . 



10 



V.P. Gerdt, Yu.A. Blinkov and V.V. Mozzhilkin 



Algorithm: HeadReduce(Q, T, -<) 



Input: Q and T, sets of triples; -<, a ranking 


Output: j7"-head reduced set Q modulo T 


1 


S:=Q 


2 


Q:=0 


3 


while 5 / do 


4 


choose p G S 


5 


S:=S\{p} 


6 


h :=HeadNormalForm(p, T) 


7 


if h ^ then 


8 


if lm(pol(p)) ^ lm(/i) then 


9 


Q := QU {h,h,0} 


10 


else 


11 


Q := Q U {p} 


12 


fi 


13 


else 


14 


if lm(pol(p)) = lm(anc(p)) then 


15 


for all € S \ anc(g) = pol(p)} do 


16 


S:=S\{q} 


17 


od 


18 


fi 


19 


fi 


20 


od 


21 


return Q 



Algorithm: HeadNormalForm(p, T, -<) 



Input: T, a set of triples; p, a triple; -<, a ranking 




Output: /i = HNFj(p,T), the jT-head normal form of pol(p) mo 


dulo T 


1 


h := pol(p) 




2 


G:={pol(g)\geT} 




3 


if lm(/t) is jT-irreducible modulo G then 




4 


return h 




5 


else 




6 


take g G T s.t. lm(/i) is ^-reducible modulo pol(g) 




7 


if lm(/i) ^ lm(anc(p)) then 




8 


if pol(p) = flo P ol(/) with f eT,9e DP(f, T) then 




9 


if Criteria(p, g) then 




10 


return 




11 


fi 




12 


fi 




13 


else 




14 


while h and 3g G T s.t. lm(/i) is ^-reducible by q := 


pol(g) do 


15 


h:=h- lc(h) 9 o ( 9 / lc( 9 )) with G T) and lm(/i) 


= 6 o lm(g) 


16 


od 




17 


fi 




18 


fi return h 
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If all the criteria are false, then the j7-head reduction of h is done by the while-loop 14-16 
in accordance with the definition of the head reduction in 1)12)) . 

The last algorithm TailNormalForm completes ^-reduction of the j7-head reduced polyno- 
mial in the input triple by performing its ^-tail reduction. This algorithm is invoked at step 13 
of the main algorithm GrobnerBasis. The tail reduction is performed in the while-loop as 
a sequence of elementary reductions (fTTj). 



Algorithm: TailNormalForm(p, T, -<) 



Input: p, a triple such that pol(p) = HNFj(p,T); T, a set of triples; -<, a ranking 


Output: h = NFj{p, T), the (full) ^-normal form of pol(p) modulo T 


1 


G:={pol(g) | g eT} 


2 


h := pol(p) 


3 


while h has a term t = ai? o y 3 which is ^/-reducible modulo G do 


4 


take g G G s.t. $ o y 3 = 9 o lm(g) 


5 


h :=h-\c{h)$o{g/\ C {g) 


6 


od 


7 


return h 



Because of the lack of an appropriate collection of benchmarks for linear finite-difference 
polynomial systems, the algorithmic efficiency of algorithm GrobnerBasis can be indirectly 
analyzed by running its polynomial (non-difference) counterpart |27|I33| for the extensive bench- 
marks collection in [381 139j . Some timings for our polynomial implementation can be found on 
the Web page |2H|. 

Recently, the algorithm in its difference version was implemented in Maple |29| . Just this 
implementation was used for generation of linear finite-difference schemes as described in the 
next sections. Though one needs special and intensive benchmarking for linear difference sys- 
tems, our first experimenting with the Maple implementation and with that for commutative 
polynomials gives us a good reason to expect that the following merits revealed for the pure 
polynomial version |57j hold also for the difference one: 

• automatic avoidance of some useless reductions; 

• weakened role of the criteria: even without applying any criteria the algorithm is reaso- 
nably fast; 

• smoothed growth of intermediate coefficients; 

• fast search of a reductor which provides the elementary Janet-like reduction IJ11JI of a given 
term. It should be noted that there can be at most one reductor |27j : 

• natural and effective parallelism. 

5 Illustrative examples of PDEs 
5.1 Laplace equation 

In this section we illustrate the approach of Section 2 to the automatic generation of difference 
schemes by simplest elliptic, parabolic and hyperbolic equations. To compute Grobner bases 
providing the elimination of the partial derivatives to construct difference schemes we used the 
Maple package |29| implementing the algorithms described in the previous section. 
We start with the Laplace equation 13 El El El 



(17) 
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and rewrite it as the conservation law l[T)) 

—Uydx + u x dy = 0. (18) 



r 



Now we add the relations (JSJ) for the partial derivatives u x and 

CDk+2 



u x dx = u(x j+2 ,y) - u(xj,y), / u y dy = u(x,y k+2 ) - u(x,y k ). (19) 



.'//,■ 



Thus, we obtain the system of three integral relations (|18|) . (|19j) for three functions 

u(x,y), u x (x,y), u y (x,y). 

To discretize this system we choose the rectangular contour of Fig. ^ on the orthogonal and 
uniform grid (J2J) with 

h 1 = h 2 = h (20) 
and use the midpoint integration method for both 1)18(1 and (|19|l . This yields the system: 

-((Uy) j + lk ~ (Uy) j + lk + 2 ) + {{U x ) j + 2k + l ~ (Uy) jk+ l) = 0, 

{u x )j + lk ■ 2/l = Uj +2k - Uj k , 

( u y)j fc+i • 2 /i = Uj fe +2 — fc- 
Rewritten in terms of difference polynomials in the ring Q){u,u x ,u y } (see Section 2) it reads: 

{0 X e 2 y - X ) O Uy + (6 2 x 6y - 6y) O = 0, 

2y i o Uj .-(^-l)ou = 0, 
o Uy - (0 2 , - 1) ou = 0. 

Computation of the Grobner basis for the elimination ranking (Section 2) with u x y u y y u and 
6 X y 6 y gives: 



6 X o Ux - -^(6 2 x - 1) o u = 
1 

,y ^ U, x -T V X ^ U,y — ~~ X ; "T V^J/ 



o Ux + 9 X o Uy - — {9 x 6 y ((6 2 x - 1) + {6 2 y - 1))) o « = 0, 



01 ° u y - ^ (^((^ - 1) + (0j - 1)) - fl# a 2 - 1)) o u = 0, 
u y - jj^(0 2 - I) o u = 0, 

i (<£flj + ^fl} _ 40^2 + ^2 + 02 } Q „ = Q 

The latter equation with eliminated u x and it,, is the standard difference scheme with the 
central approximation of the second-order derivatives in (|17j) written in double nodes 

u j+2k — %Uj_k + u j-2k Uj k+ 2 — 2Uj k + Uj k ^2 _ „ 

Similarly, the trapezoidal integration rule for relations ()19|) generates the same difference 
scheme but written in ordinary nodes 

u j+l k — 2Uj k + fc ttj fc + i — 2Uj k + Uj fc_j _ 
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n+l 




3 3 + 1 3 + 2 

Figure 2. Integration contour for heat equation. 



5.2 Heat equation 

Consider now the heat equation IU El El HI IE] 

u t + au ra = 
in its conservation law form 

-au x dt + udx = 0. 



(21) 



The integrand in (|21jl contains the only partial derivative u x . Hence we add the single integral 
relation 



Xj + 1 



u x dx = u(xj + i, t) — u(xj, t). 



(22) 



Again we discretize u(x,t) and u x (x,t) on the orthogonal and uniform grid with the spatial 
mesh step h and the temporal mesh step r, and choose the contour shown in Fig. |2 Then, 
applying the midpoint rule for the contour integral and the trapezoidal rule for the relation 
integral we find two difference equations for two indeterminates u, u x in the form 

a- (i + e t -9 2 x - e t e 2 x ) OUx -2h (e x e t -e x )o U = o, 

~ (9 X + l)ou x -(6 x -l)ou = 0. 

By elimination of u x by means of the Grobner basis with u x >- u we obtain the famous 
Crank-Nicholson scheme H © □ El El 



0. 



r 2/i 2 
The same scheme is obtained for the midpoint integration method applied to 1)22(1 . 

5.3 Wave equation 

The wave equation [3 il El El El El 

u t t ~u xx = 

in the conservation law form is given by 



u x dt + utdx = 0. 
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Choosing the same grid with the mesh steps (|2Uj) . the contour of Fig.^ an d integral relations ()19|) 
as are used in Section 5.2 for the Laplace equation and applying the midpoint rule for 
the contour integral and the trapezoidal rule for the integral relations we obtain the operator 
equations 

(e x - e x e 2 t ) o Ut + (6% - o t ) o Ux = o, 

~ (9 X + l)ou x -(6 x -l)ou = 0, 

^ (0 t + 1) o ut - (p t - 1) o u = 0. 
The Grobner basis method yields the standard difference scheme 
u] +1 + u r ;~ l - U ] +1 - = 0. 

5.4 Advection equation 

Consider now a simple form of the Advection (or convection or one-way wave) equation [HI 13 

Ut + vu x = 0, v = const. (23) 

Being of first order, the equation Q23|) has already the conservation law form By this 
reason, to generate a difference scheme we shall not convert the equation into the integral 
form (J2J). In the latter case one has nothing to eliminate. Instead, we consider equation (|23j) 
together with the integral relations: 

u t + vu x = 0, 

r-t-2 rX2 

Utdt = u(t2, x) — u(t±, x), / u x dx = u(t, X2) — u(t, x\). (24) 

H\ J x\ 

Discretization of u, ut and u x on the orthogonal and uniform grid with the mesh steps h and r 
and the explicit integration formula for the upper integral relation in (|24ft together with the 
midpoint integration rule for the lower relation give the difference system: 

ut + v u x = 0, r ut = (9t — 1) o u, 2h9 x o u x = (9 2 — 1) o u. (25) 

Let us apply the operator 9 X to both sides of the middle equation in (|25[) and then use the 
Lax method, that is, replace 9 X with (9 X + l)/2 in the second term of the right-hand side. This 
replacement yields 

u t + vu x = 0, 

0x°UfT-( t x - J U = 0, 

20s ou x ■ h - {9 2 x - 1) o u = 0. 
The lexicographical Grobner basis for the elimination ranking ut >- u x y u with 9t >- 9 X , is 
u t + vu x = 0, 

2h0 x o Ux - (6l - 1) on = 0, 

(2h 9 x 9 t -h{9 2 x + l) + t{9 2 x -1)-u)o U = 0. 

Its last element gives the scheme: 



U 



"2 2h 
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6 Burgers equation 

6.1 Conservation law form 

Consider Burgers equation |H EJ |S1 E] in the form 

u t + fx = vuxx, v = const, (26) 

where we replaced u 2 by the flux function / in order to avoid computation of nonlinear difference 
Grobner bases, v is called the viscosity. This equation exhibits some difficult features from the 
point of view of simple finite difference schemes due to the term / = u 2 . Let us, first, convert 
equation (|26|) into the conservation law form 

(vu x — f)dt + udx = 0. 

Then choose the contour of Fig. ^ and add the integral relation 

u x dx = u(Xj + 2, t) — u(xj, t). 

~3 

Denoting as above the temporal and spatial mesh steps by r and h, and applying the midpoint 
integration rule we obtain the system: 

h (M? -e x )o U -T (e% - e t ) o (uu x - f) = o, 

2h9 x ou x - (6l - 1) on = 0. 
Its Grobner basis form for the elimination order with u x y u y f and 0t >~ e x is given by 

2 urh e t o Ux + 2 h 2 e x (ef -i)o«+2 rhe t (el - 1) o / - u T e t e x (el -i)o U = o, 
2he x o Ux - (el - 1) ou = o, 

2 h 2 e 2 x (e 2 t -\)o U - vre t (e x -2^ + i)o M +2r/i e t e x (el - 1) o / = o. 

The obtained difference scheme 

n+2 _ n n+l n„.n+l , n+1 rn+1 f n+l 

u j+2 u j+2 _ ^U j+2 +U j J j+3 J j+1 _ 

t 2h 2 + h ~ [ ' 

is the standard explicit scheme with forward time and forward space differencing. It is well- 
known that schemes of this type are unstable [5JIB1IH]- Furthermore, by using implicit schemes 
one can provide the von Neumann stability 1 . However, all such schemes are usually not very 
satisfactory when one considers non-smooth or discontinuous solutions (shock waves) of Burgers 
equation. 

6.2 Lax method 

To exploit more flexibility and freedom in our difference elimination approach to generation of 
finite-difference schemes, we go back to the original differential equation (|26|) and consider it 
together with the integral relations providing the elimination. For discretization of the relations 
we combine the midpoint rule for integration over x with the explicit integration over t and 
apply the Lax method to the last integration: 

u t + f x = vu xx , (utYj + (fx)] = v(u xx ) n j , 

x For example, if one uses the central differencing in the last term of l|27|l . and the Crank-Nicolson approach 
to the second (diffusion) term. 
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U t dt = U, U t T = U^- Ul+2 ^ U \ 

kdx = f, 2 h(f x )? +1 = f? +2 - (28) 

u x dx = u, 2 h(u x )] +1 = u] +2 - uj, 

u xx dx = u x , 2h{u xx Yj + i = ~ { u x) 1 j ■ 

The Grobner basis based elimination with u xx y ut y u x y f x y u y f yields the scheme 

2~r + 2h " IP °" (29) 

One can also use the trapezoidal rule for the spatial integrations. This derives other schemes. 
Since there are three spatial integrals in (|28|). by selecting either the midpoint or the trape- 
zoidal rule for these integrals, we obtain eight possible variants of the difference schemes. Our 
computation with the Grobner bases reveals seven different schemes. Apart from (|29j) there are 

2(n^ 2 1 + ^+ 1 )-(n- +3 + ^ +2 + ^ +1 + ^) (/;•., • /;'.,) (/;•., • /;■) 

4r Ah 

l , K +3 -^ +2 )-K + i-^) m) 



2t 2h h 2 

2{u% + 2<+ 1 + u]+l) - (u? +4 + 2u] +3 + 2u] +2 + 2u» +1 + 



8r 

+ 8h " P ' {6 ] 



2(^+ 1 +^+ 1 )-(n J " +4 + ^ + 3 + ^ +2 + ^ +1 ) ,/;;-. :{ /;'., 

4r h 

((u] +5 + u] +4 ) - 2u- +3 ) - (2u] +2 - (u] +1 + up) 



8/i 2 



(33) 



(34) 



2(n^+u^)-(^ +3 + ^ +2 + ^ +1 + ^) ff +2 -ff +1 

4r /i 

K + 3-^ +2 )-(n? +1 -n?) 

" 2/> 2 

2(tff 3 1 + 2u n +l + iff* 1 ) - (n? +4 + 2^ +3 + 2u] +2 + 2n? +1 + «?) - 

8r 2/i 
= /i + 3-2^- 2 + ^ +1 (35) 

Just the scheme (|34|) is obtained twice in the course of generating eight schemes. 



6.3 Two-step Lax— Wendroff method 

Our Grobner basis based technique can also be applied to generate other types of difference 
schemes. For example, one can generate two-step Lax- Wendroff schemes |41| . Let u and / 
denote the values of u and / on the intermediate time levels. Then, applying again the midpoint 
rule for the spatial integrals, gives the following difference system: 

Utj+fxj = V Uxx], 
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„+l U ]+2 + ul j 



U tj T = Uj 

2fx] +1 h = f? +2 -f?, 

% u x j+l h = Uj+2 ~ u j i 

21*^2; h = U x j_(_2 U x j ; 
— n 1 ~~t n . . — n 



U U + fxj = v U 



2 fx j+l h = fj+2 ~ fj 

2u x ] +1 h = u] +2 -u 



■n 



1u x X h — 1l x j_(_2 U x j . 

For the elimination ranking with 

u xx yu xx y u x yn x >- u t yut y f x >- f x >- f >- u >- f ^ ^ 

the Grobner basis contains the Lax-Wendroff scheme 

U ]+2 ~ ( n j+3 + U ]+l) /J+3 ~ /J+l = ~ 2u i+2 + U i 



2r 2/i 4/i 2 

^+3 ~ ^+2 7j+3 ~ / j+l _ gjj-4 ~^+2 +^ ,„ R s 

27 + 2^ " ^ IP ' W 

With all possible combinations of the trapezoidal and midpoint rules one obtains 49 different 
Lax-Wendroff schemes. 



6.4 Differential approximation 

To analyze properties of a difference scheme it can be useful to compute its differential ap- 
proximation [30J that is often called the modified equation(s) of the difference scheme. There 
are whole classes of different schemes for which their stability properties can be obtained with 
the aid of the differential approximation [2]. For all that, in many cases, the computation can 
be easily done with modern computer algebra software. In our research we use Maple |17j . 
Consider, for example, the schemes (|29*|l - (j55)l . 

Their differential approximation for f = u 2 and with collection of the coefficients at r, h 2 , 
h 2 /t is given by: 

+ (*) h +t; u xx— + ■■■ ■ 

Z T 

The schemes (|2*9*|) - (|35J) differ in the coefficient (*) at h 2 only. These coefficients are as follows 
1 1 1 

("j2J — U XXXX V ~\~ 77 V'xxx'U* ~i~ 77 ^ra^ii 

6 2 

1 1 1 

77 ^xxxx^ ~H T77 U X xxU 77 ^xxU x , 

12 4 

5 1 1 

("j4J — u xxxx v + — u xxx u + — u xx it X: 

1 1 

("J"J 7^ "rara^ 7^ V'XxxU U XX U X , 

D 
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V ) 12 xxxx 3 xxx 2 XX X, 

1 1 1 

J 77 Uxxxx^ ~\~ 77 'UxxxU ~\~ 77 

2 

1 1 1 

("JoJ — Uxxxx^ ~T V'XxUx 777 U XXX U. 

o 4 12 

Thereby, comparison of differential approximations for schemes ()29|) - (|;-i5|) shows that 

• all the schemes provide the same order of approximation in r, h; 

• they have identical linear numerical dissipation (viscosity) [HE] determined by u xx h 2 / (2 r); 

• the schemes possess similar dispersion properties with distinction in the rational coefficients 
of the differential polynomial in u multiplied by h 2 . 

As to scheme (|27|l. the right-hand side of its differential approximation reads 

--V U XXXX + (UxxxU + 2u xx U x )v - U X U - -U U xx j T 

(I 1 1 \ 2 

~(~ I g^xxxa;^ g^xxx^ 2^ X3; ^ X ) ~\~ ' ' ' ■ 

This explicitly shows instability of the scheme which does not yield linear numerical viscosity. 

We obtained also analogous results on stability and on close properties for the different 
Lax-Wendroff schemes of type (|36|) and its variations due to the choice of different numerical 
integration rules for the spatial integrals. 

6.5 Godunov method 

It is especially difficult to simulate numerically nonsmooth and discontinuous solutions which 
are among most interesting problems in computational fluid and gas dynamics ^ HJ 03 |Sj. 
Most of the known difference schemes fail to handle these singularities. The most appropriate 
numerical approach to such problems was developed by Godunov ^ H2] and based on solving 
a local Riemann problem [3J as a cornerstone of the Godunov's scheme generation. There 
are special numerical Riemann solvers, for example |43| . designed for these purposes and for 
application to computational fluid dynamics. 

Instead of the use of numerical Riemann solvers, we apply the Grobner bases technique to 
generate the Godunov- type scheme for inviscid Burgers equation when v = in (|26|) . For this 
purpose we discretize the corresponding system in ((2*5|) in the following way 

u t ] + fx] = 0, 
u t ]r = u] +1 -u], 

(fx?h- (f* +1 - f?))(f m * +1 h- (/JVi " #)) = o, 
2 u x " +1 h = itjl|_2 — u j i 

2 u xx h = u x ™_|_2 — u x j • (37) 

Here, the third equation contains in its left-hand side the product of two different solutions 
for the flux function / of the local Riemann problem [33]. Therefore, we add to the system 
composed of the original differential equation and discrete forms of the integral relations for 
partial derivatives Ut, u x , u xx the nonlinear difference equation on / and f x containing solutions 
of the local Riemann problem. 

Since the Riemann condition on the flux is now a constituent of the difference system, an 
elimination of all the partial derivatives of u and / gives the difference scheme consistent with 
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that condition. To do the elimination from the nonlinear system (|37j) we apply the Grobner 
factoring approach [33]: if a Grobner basis contains a polynomial which factors, then the com- 
putation is split into the computation of two or more Grobner bases corresponding to the factors. 
In doing so, we choose the elimination ranking 

u xx y u x y m y fx y f y u 

and compute two Grobner bases, for every factor in (|37|). Then we compose the product of two 
obtained difference polynomials in u and / that gives us the Godunov-type difference scheme: 

,,H+1 _ „.n fn fn \ /,, n +l _ ,,n fn fn \ 

U j+2 U j+2 Jj+2 Jj+1 \ I U j+2 U j+2 Jj+3 Jj+2 \ , . 
T + h ) • ^ T + h J' ' (38) 

Below (Section 7) we compare schemes (|'29|) . ()36(l and (|38)1 by numerical simulation of a dis- 
continuous solution. 



7 Falkowich— Karman equation 

Consider now the nonlinear two-dimensional Falkowich-Karman equation |3JJ] 

<Pxx(K - (7 + l)<Px) + Vyy = (39) 
describing transonic flow in gas dynamics in its non-stationary form 

<p xx (K - (7 + l)<Px) + Vyy ~ tyxt ~ Vtt = 0. (40) 

This form can be used to find a stationary solution by the steady-state method. We rewrite 
equation (|4U|) into the conservation law form 



[tn+l / r / ( 7 + 1) \ \ fXj+2 fVk+2 

J \^J> -Vydx + ipx \ K — ip x J dyjdt-J J (2ip x + ipt) 



tn + l 

dxdy = 0. 

t n 



Here we use a grid decomposition of the three-dimensional domain in (x, y, t) into elementary 
volumes. Fig. |21 shows an elementary volume. 

Again we add the integral relations for partial derivatives with the use of the trapezoidal 
integration rule for (p x , (p y and the midpoint rule for (ft. 

Then we obtain the nonlinear operator equations: 

ie x -e x el)oy y + (ele y -o y )oL x (k-^^-^A yj - 2hr 

- (0 t - l)(6 2 x e y -9 y )o2ip-2h-e x e y oip t .4:h 2 = 0, 

(0 X + 1) O Lf x ■ ^ = (0 X -l)o if, 
(0 v + l)o(py~ = (0 y -l)o if, 

9 t o Vt .2T = (9?-l)o V . (41) 

Because of nonlinearity in the initial differential equation 1)40(1 , the difference system obtained 
is also nonlinear. By this reason the Maple package [21] implementing algorithm GrobnerBasis 
(Section 4) is inapplicable to (|41jl. Since there is no software available for computing difference 
Grobner bases for nonlinear systems, we had to perform hand computations in accordance with 
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n + 1 




k + 2 



k + 1 



J + 2 k 

Figure 3. Cell for Falkowich-Karman equation. 



the above described algorithm and with an assistance of Maple to check some intermediate 
results. 

In these calculations we used the lexicographical ranking such that tp x y <p y y ip% y <p and 
8 X y y y Of The resulting Grobner basis has the form: 



(B x - l) 2 y o tp ■ (7 + l)Vxr + QxOy °<Pfh 



((B x - lfe y o<p.(K 



ife y 0(7 + iv) + e x {e y -ifo <p)- 



] v - e tt )(e t - 1) o <p, 



(By + 1) O <p y = (By - 1) O <f ■ ~, 



h ° ft 



■ K 



&l + B x - \)B y B t o 



B X (B X - l) 2 B y B t o <p ■ [((^ - l) 2 ByB t o <p ■ ( 
+ B x (B y - l) 2 8 t o^l-(B 2 x - l)B y {B 2 t -0 t )o<p- 8 x B y (B t - l) 2 o <p • A 



(7 + 1) 



+ (B x - lf&yB t o ip ■ [B X (B X - l) 2 ByB t o cp ■ ( K 



f x + B x -l)B y B t o^±}l ip 



i\2 Q „ ,A T ro3 a\a(a1 a \ n ,„ c&aia i\2 



+ 8i(B y - iy& t o <p) - - (8 X - & x )8 y (&t -e t )o<p- BiBy{B t - 1)' a <p 



The last element is the finite-difference scheme for equation (|40j) : 



+ <^_ ljt - ^_ 2fc ) + {Vj-lk+1 ~ 2 <Pj-lk + <Pj-lk-lj) 



n-1 \ h 



+ {^ k - 2^_ lk + vUk) ■ [((¥>?+!* " + <fi-ik)(K ~ ^T^l+lk ~ 



0. 



(42) 



By construction, this scheme is fully conservative and does not involve switches that is typical 
for computing transonic flows [4*7| , 
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fc + 1 



k - 1 

3-2 j I j J+l 

Figure 4. Stencil for stationary Falkowich-Karman equation. 



In its stationary form scheme (J12 
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0. 



(43) 



is related to equation (J3U). Here symbols D x and are the forward differencing operators 
and D xx and D yy are the central second-order differencing operators with respect to x and y. 
The stencil for scheme (j43j) is shown in Fig. 

It should be noted that, unlike the original differential equation (|40j) which is quadratically 
nonlinear, both schemes ()42|) and l)43|) have the the cubic nonlinearity in the grid function. This 
is the rigorous algebraic consequence of the difference system ()41|) . In accordance to the well- 
known fact [inj; that algebraic elimination of variables from a nonlinear system leads generally 
to increase of its degree of nonlinearity. 

As an application of this scheme, in the next section we consider an example of one-dimen- 
sional transonic flow with shock- wave taken from |47| . 



8 Numerical experiments 
8.1 Burgers equation 

We used schemes (|29j) . (|36[) and (j38j) for numerical simulation in < x < 1 of the following 
Riemann problem for the inviscid Burgers equation (|26|) with v = 



u t + u x u = 0, 
and discontinuous initial condition u(x,0): 



u(x, 0) 



Ui, 



< x < 1, 

1 < x < 1. 



(44) 



(45) 



In ()45j) the initial data at t = is a piecewise-constant function with the state ti/ on the left of 
the discontinuity x = and the state u r on the right of the discontinuity. We consider v = 0, 
since in this case the problem (j44j) . ()45|) admits the exact solution: 



' 1 u z + n r 
u(x,t) = u t H ( - H — t 



I + u r H [ x l — — - 1 



(46) 
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I 1 1 | 1 1 

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

Figure 5. Lax scheme (12911 with Courant num- Figure 6. Lax scheme (|29|l with Courant num- 
ber 0.9. bcr 0.1. 




Figure 7. Lax-Wendroff scheme (j86(l with Cou- Figure 8. Lax- Wcndroff scheme (|36J) with Cou- 
rant number 0.9. rant number 0.1. 



Here H (y) is the Heaviside step function |48| whose derivative is the Dirac delta function: 
r y<0, 
H(y) = l \ y = 0, ± H (y) = 8(y). 

{ i y>o, 

Physically, the solution ()46j) defined by the initial condition (|45|) represents a shock wave which 
moves with constant speed (ui + u r )/2 without changing its shape. 

In our numerical simulation the values of u\ and U[ were chosen as 0.8 and 0.2. The pictures 
below demonstrate the numerical solution of the Riemann problem ()44|) . ()45|) at time t = 2/3. 
Solid line shows the exact solution ()46|) . and the numerical results are depicted by green dots. 
For the ratio r/h of mesh steps which is called Courant (or Courant -Friedrichs-Levy) number 
we have chosen the two values 0.9 and 0.1. 

All schemes are numerically stable. For schemes (|29|) and (|3(i[) their stability is analytically 
showed by the differential approximation (Section 6.4). Because of the nonlinearity in f x in the 
third equation of Godunov scheme (|38|) we did not compute the differential approximation for 
this scheme. 

As can be expected (cf. (31), the dispersion effects in schemes (|29|) - (|36[) become stronger for 
the smaller value of the Courant number (Figs. E] and |HJ). Qualitatively [491 150j . the behavior 
of Lax scheme (|29j) in Figs. El is typical for the classical first-order schemes when they are 
applied to problem (|44p - ()45|) whereas the Lax-Wendroff scheme (|36|l behaves as the second- 
oder method. The Godunov scheme ()38|) . as a shock capturing one (cf. |49U50| L is much better 
for numerical description of solution (|4ri|) than the schemes (j29j) - (|3(ij) . and does not reveal its 
sensitivity to the value of the Courant number. 
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Figure 9. Godunov scheme (|38fl with Courant Figure 10. Godunov scheme l|38fl with Courant 
number 0.9. number 0.1. 





Figure 11. Initial numerical approximation for 
equation Ij39(l . 

8.2 Falkowich— Karman equation 



Figure 12. Numerical solution of equation l|39|) . 



Now we consider the application of difference scheme (|43[) to the one-dimensional stationary 
transonic flow in a channel with a straight density jump |47j . The exact shock- wave solution 
of equation 1)39(1 at < x < 1 is shown in Figs. ^2 and ^] by solid red line. Circles depict the 
numerical data obtained from difference scheme (|43|) . As an initial approximation, the parabola 
was chosen satisfying the following boundary conditions: at the left, both the function and its 
derivative are fixed by the values from the exact solution; at the right, the only function is 
bound to the exact solution. 

As one can see from Fig. 1121 scheme (|43|l possesses a stable and uniform convergence to the 
exact shock-wave solution. Because, by its construction, the scheme is fully conservative, it does 
not reveal non-uniqueness of solutions that is typical for the traditional difference schemes |47| . 

Moreover, the size of the shock transition zone is just one spatial mesh step that is a conse- 
quence of preserving at the discrete level of all algebraic properties of the initial PDE (|4Uj) . This 
is a result of algebraic difference elimination provided by the Grobner bases method. Another 
merit of scheme (|43|) is that it does not involve switches that is typical for computing transonic 
flow as we already pointed out in Section 7. 

This example shows a principal possibility of constructing difference schemes for transonic 
flow without switches and with the same stencil for both subsonic and supersonic flow. 



9 Conclusion 



In the present paper we have shown that the Grobner bases method, being a universal algorithmic 
tool for linear difference algebra, can be effectively applied to the construction of differences 
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schemes for linear PDEs with two independent variables and with rational function coefficients. 
Owing to the Grobner bases, this construction is an algorithmic procedure. It consists in elimi- 
nation of partial derivatives from the system of difference equations composed from a discrete 
version of the original PDEs (on an orthogonal uniform grid) and numerically approximated 
integral relations between the unknown functions and their partial derivatives. As this takes 
place, the difference scheme obtained may depend on the choice of the integration contour and 
numerical approximations for integral relations. 

The method is especially efficient when a PDE or a system of PDEs admits the conservation 
law form. In this case the difference schemes obtained are fully conservative. The structure of 
a scheme generated may depend on the choice of integration contour and numerical integration 
rules. In so doing, it is not clear a priori which integration rule leads to a better scheme. 

We also described an efficient algorithm for the construction of Grobner bases for linear 
difference ideals. The algorithm is based on the concept of Janet-like reductions. Its first 
implementation in Maple is already available, and we used this implementation in the generation 
of all linear difference schemes presented in the paper. 

For classical linear PDEs such as the Laplace equation, the Heat equation, the Wave equation 
and the Advection equation our algorithmic technique leads to the well-known finite difference 
schemes. For Burgers equation we generated several schemes based on the Lax and Lax- 
Wendroff methods and computed their numerical dissipation and dispersion by the differential 
approximation (modified equation) method. By example of Burgers equation we also demon- 
strate that it is possible to combine the Godunov method with Grobner bases to derive a shock 
capturing scheme. 

The non-traditional cubic nonlinear difference scheme generated by our difference elimination 
method for the Falkowich-Karman equation describing transonic flow in gas dynamics possesses 
a number of attractive properties in comparison with traditional schemes. Among them there 
are a stable convergence in time to the exact solution with a one-dimensional shock wave and 
absence of switches. It should be noted, however, that due to its cubic nonlinearity, scheme Q43|) 
has a slower convergence in comparison with the traditional schemes specially optimized for 
numerical simulation of transonic flows in gas dynamics. By this reason one needs additional 
research for optimizing nonlinear schemes obtained by the difference elimination. 

As we already mentioned in the introduction, algorithm GrobnerBasis admits a generali- 
zation to polynomial-nonlinear systems of difference equations exactly in the same way as the 
differential involutive algorithm of paper (221- In doing so, if every equation in the initial system 
is linear with respect to the highest ranking difference term and this property of the system 
is not violated during its completion to involution, then algorithm GrobnerBasis will work 
correctly and provide the desirable output. Such is indeed the case for system (|4"Tj) . In the most 
general case of a difference system with polynomial nonlinearity, it can be split into a finite 
number of subsystems such that every subsystem can be converted into the Grobner basis form 
by applying our algorithm. The underlying splitting algorithm is a difference analogue of that 
described in |31| . The latter algorithm is similar to the splitting algorithm implemented in the 
library package diffalg in Maple. 

The above described approach can be also generalized to PDEs with three and more indepen- 
dent variables. Thus, if PDEs admit the conservation law form, then one can use multidimen- 
sional analogues of equations Q and (J2J) together with their elementary volume discretization. 
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